﻿//Mexican Hat wavelet
namespace Diplom.Wavelets
{

    using System;
    using System.Collections.Generic;
    using System.Linq;
    using System.Text;
    using Diplom.NumericArrays;
    using System.Drawing;

    public static class Morlet
    { 
        private static Double Base = 2.0;
        private static Double Move = 1;
        private static Double Norma = Math.PI/Math.Log(Base,Math.E);        

        /// <summary>
        /// Wavelet function
        /// </summary>
        public static double MorletWavelet(double x)
        {
            //double res = (1-x*x);//MHAT
            //res *= Math.Exp(-x*x*0.5);

            double res = Math.Exp(-x * x / 2.0) * Math.Cos(2 * Math.PI * x); //Real part Morlet

            //double res = -x * Math.Exp(-x * x * 0.5); //WAVE

            //double res = Math.Exp(-x * x * 0.5) - 0.5 * Math.Exp(-x * x / 8.0); //DOG

            /*if (x < 0) return 0; //HAAR
            if (x < 0.5) return 1.0;
            if (x < 1) return -1.0;
            return 0.0;*/

            /*if (Math.Abs(x) < 1e-6) return 0.0; //LP-LittleWood_Paley
            else return(Math.Sin(2*Math.PI*x)-Math.Sin(Math.PI*x))/(Math.PI*x);*/

            /*if (x < -1.0) return 0; //FHAT
            if (3*x<-1.0) return -0.5;
            if (3 * x < 1.0) return 1.0;
            if (x < 1.0) return -0.5;
            return 0.0;*/

            return res;
        }

        /// <summary>
        /// Wavelet transform
        /// </summary>
        public static Numeric2DArray Transform(ref Numeric1DArray source)
        {
            int m = 0;
            double cur = 1;
            while (cur <= source.Length)
            {
                m++;
                cur *= Base;
            }
            cur = source.Length;
            Numeric2DArray waveCoeff = new Numeric2DArray(m + 1, (int)cur + 1);
            for (int i = 0; i <= m; i++)
            {
                double a = 1 / Math.Pow(Base, i + 0.0);
                double sqrt_a = Math.Sqrt(a);
                for (int j = 0; j <= (int)cur; j++)
                {
                    double res = 0;
                    for (int k = 0; k < source.Length; k++)
                    {
                        double phi = MorletWavelet((k - j*Move)*a);
                        res += (source.Items[k] * phi);
                    }
                    res *= sqrt_a;
                    waveCoeff.Insert(i, j, res);
                }
            }
            waveCoeff.InitCharacters();
            return waveCoeff;
        }

        /// <summary>
        /// Wavelet analys of signal in logariphim scale
        /// </summary>        
        public static Numeric2DArray Analys(ref Numeric1DArray source, int NA, int NB,bool Modulate,bool Remove)
        {
            int m = 1;
            double cur = 1;
            while (cur*Base<= source.Length)
            {
                m++;
                cur *= Base;
            }
            
            double Amin = 0;
            double Amax = m;
            double delA = (Amax / (NA - 1));
            double Bmin = 0;
            double Bmax = source.Length;            
            double delB = (Bmax-Bmin) / (NB - 1);
            Numeric2DArray waveCoeff = new Numeric2DArray(NA, NB);

            double cura = Amin;
            for (int i = 0; i < NA; i++)
            {
                double a = 1 / Math.Pow(Base, cura);
                double sqrt_a = Math.Sqrt(a);

                double b = Bmin;
                //int Ja = (int)(1.08/(delB*a));
                int Ja = (int)(1.08 *cura/ (delB));
                for (int j = 0; j < NB; j++)
                {
                    if (Remove)
                    {
                        if (j <= Ja || j >= NB - Ja)
                        {
                            waveCoeff.Insert(i, j, 0.0);
                            continue;
                        }
                    }
                    double res = 0;
                    for (int k = 0; k < source.Length; k++)
                    {
                        double phi = MorletWavelet((k - b) * a);
                        res += (source.Items[k] * phi);
                    }
                    res *= sqrt_a;
                    if (Modulate)
                    {
                        res = Math.Abs(res);
                    }
                    waveCoeff.Insert(i, j, res);

                    b += delB;
                }
                cura += delA;
            }
            waveCoeff.MinTime = Bmin;
            waveCoeff.MaxTime = Bmax;
            waveCoeff.MinFreak = 0;
            waveCoeff.MaxFreak = m;            
            waveCoeff.InitCharacters();
            return waveCoeff;
        }

        /// <summary>
        /// Wavelet analys of signal in logariphim scale
        /// </summary>        
        public static Numeric2DArray Analys(ref Numeric1DArray source, int NA, int NB, bool Modulate, bool Remove,RectangleF analizeRect)
        {
            #warning это перегруженный вариант вейвлета, для заданного кусочка, его надо испоьзовать отдельно
            #warning дополнительный параметр - RectangleF (X - начало прямоугольника по времени, Width - ширина временного прямоугольника, Y - начальная частота, Height - высота частотного прямоугольника)   

            double Amin = Math.Log(analizeRect.Y, 2.0);
            double Amax = Math.Log(analizeRect.Y + analizeRect.Height, 2.0);
            double delA = ((Amax-Amin) / (NA - 1));
            double Bmin = analizeRect.X;
            double Bmax = analizeRect.X + analizeRect.Width;
            double delB = (Bmax - Bmin) / (NB - 1);
            Numeric2DArray waveCoeff = new Numeric2DArray(NA, NB);

            double cura = Amin;
            for (int i = 0; i < NA; i++)
            {
                double a = 1 / Math.Pow(Base, cura);
                double sqrt_a = Math.Sqrt(a);

                double b = Bmin;
                //int Ja = (int)(1.08/(delB*a));
                int Ja = (int)(1.08 * cura / (delB));
                for (int j = 0; j < NB; j++)
                {
                    if (Remove)
                    {
                        if (j <= Ja || j >= NB - Ja)
                        {
                            waveCoeff.Insert(i, j, 0.0);
                            continue;
                        }
                    }
                    double res = 0;
                    for (int k = 0; k < source.Length; k++)
                    {
                        double phi = MorletWavelet((k - b) * a);
                        res += (source.Items[k] * phi);
                    }
                    res *= sqrt_a;
                    if (Modulate)
                    {
                        res = Math.Abs(res);
                    }
                    waveCoeff.Insert(i, j, res);

                    b += delB;
                }
                cura += delA;
            }
            waveCoeff.MinTime = Bmin;
            waveCoeff.MaxTime = Bmax;
            waveCoeff.MinFreak = Amin;
            waveCoeff.MaxFreak = Amax;
            waveCoeff.InitCharacters();
            return waveCoeff;
        }


        /// <summary>
        /// Wavelet syntes
        /// </summary>        
        public static Numeric1DArray Syntes(ref Numeric2DArray source, int size)
        {
            Numeric1DArray signal = new Numeric1DArray(size);
            
            for (int k = 0; k < size; k++)
            {
                double res = 0;
                double olda = 1.0;
                for (int i = 0; i < source.NSize; i++)
                {
                    double a = 1.0 / Math.Pow(Base, i + 0.0);
                    double sqrt_a = Math.Sqrt(a);
                    double current = 0;
                    for (int j = 0; j < source.MSize; j++)
                    {
                        double phi = MorletWavelet((k - j * Move) * a);
                        current += (source.Items[i][j] * phi) * sqrt_a;
                    }
                    res += (current * a * a);// *(1.0 / a - olda);
                    olda = 1.0/a;
                }
                res /= (Norma);
                signal.Insert(k, res);
            }
            signal.InitCharacters();
            return signal;
        }

    }
}
